suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(ChIPseeker))
suppressPackageStartupMessages(library(clusterProfiler))
suppressPackageStartupMessages(library(argparse))

parser <- ArgumentParser()
parser$add_argument("-i", "--peak_file", help="peak file")
parser$add_argument("-o", "--anno_file", help="annotate peak file")
parser$add_argument("-s", "--species", default="hs", help="species hs|mm|rn (default='hs')")
args <- parser$parse_args()

if (args$species == "hs") {
  suppressPackageStartupMessages(library(TxDb.Hsapiens.UCSC.hg38.knownGene))
  txdb <- TxDb.Hsapiens.UCSC.hg38.knownGene
  annoDb <- "org.Hs.eg.db"
} else if (args$species == "mm") {
  suppressPackageStartupMessages(library(TxDb.Mmusculus.UCSC.mm10.knownGene))
  txdb <- TxDb.Mmusculus.UCSC.mm10.knownGene
  annoDb <- "org.Mm.eg.db"
} else if (args$species == "rn") {
  suppressPackageStartupMessages(library(TxDb.Rnorvegicus.UCSC.rn6.refGene))
  txdb <- TxDb.Rnorvegicus.UCSC.rn6.refGene
  annoDb <- "org.Rn.eg.db"
}

peakAnno <- annotatePeak(args$peak_file, tssRegion=c(-3000, 3000), TxDb=txdb, annoDb=annoDb)

peakAnnoList <- lapply(files, annotatePeak, TxDb=txdb, tssRegion=c(-3000, 3000), verbose=FALSE)

plotAnnoBar(peakAnnoList)

plotDistToTSS(peakAnnoList)

genes = lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
names(genes) = sub("_", "\n", names(genes))
compMF <- compareCluster(geneCluster   = genes, 
                         fun           = "enrichGO",
                         ont           = "MF",
                         pvalueCutoff  = 0.05, 
                         pAdjustMethod = "BH")
plot(compMF, showCategory = 20, title = "Molecular Function Enrichment")

genes= lapply(peakAnnoList, function(i) as.data.frame(i)$geneId)
vennplot(genes)
